home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / examples / estimate.r < prev    next >
Encoding:
Text File  |  1994-09-20  |  3.6 KB  |  175 lines  |  [TEXT/ttxt]
  1. // -------------------------------------------------------------------
  2. //
  3. // Beginning of codes to test memory allocation/deallocation.
  4. //
  5. // The following lines are in a file called `estimate.r'
  6.  
  7. mysum = function(X) 
  8. {
  9.   //
  10.   // Treat vectors, scalars and matrices differently.  If X is a scalar,
  11.   // do nothing...
  12.  
  13.   if (X.n == 1) { return X; }
  14.  
  15.   //
  16.   // if X is a matrix, return the vector of row-sums, except in the special
  17.   // case that X.nc==1 or X.rc==1; in that case, return the scalar sum.
  18.   //
  19.  
  20.   if (X.n  > 1)
  21.   {
  22.     if (X.nc == 1) 
  23.     {
  24.       tmp = 0;
  25.       for(i in 1:X.nr) 
  26.       {
  27.         tmp = tmp + X[i;1];
  28.       }
  29.       return tmp;
  30.     }
  31.     if (X.nr == 1) 
  32.     {
  33.       tmp = 0;
  34.       for(i in 1:X.nc) 
  35.       {
  36.         tmp = tmp + X[1;i];
  37.       }
  38.       return tmp;
  39.     }
  40.     tmp = zeros(X.nr);
  41.     for(i in 1:X.nr) 
  42.     {
  43.       for(j in 1:X.nc) 
  44.       {
  45.         tmp[i] = tmp[i] + X[i;j];
  46.       }
  47.     }
  48.     return tmp;
  49.   }
  50. }
  51.  
  52. initstate = function(v) 
  53. {
  54.   local (prob, tmp, i);
  55.  
  56.   tmp = 0;
  57.   prob = rand();
  58.   for(i in 1:v.nr) 
  59.   {
  60.     tmp = tmp + v[i;1];
  61.     if (tmp >= prob) 
  62.     {
  63.       break;
  64.     }
  65.   }
  66.   if (v[i;1] == 0) 
  67.   {
  68.     error("initstate: Attempting to initialize to invalid state.");
  69.   }
  70.   return i;
  71. }
  72.  
  73. nextstate = function(thisstate, trans) 
  74. {
  75.   local (prob, tmp, i);
  76.  
  77.   tmp = 0;
  78.   prob = rand();
  79.   for(i in 1:trans.nc) 
  80.   {
  81.     tmp = tmp + trans[thisstate;i];
  82.     if (tmp >= prob) 
  83.     {
  84.       break;
  85.     }
  86.   }
  87.   if (trans[thisstate;i] == 0) 
  88.   {
  89.     error("nextstate: Attempting transition to invalid state.");
  90.   }
  91.   return i;
  92. }
  93.  
  94. randwalk = function(A, P, f, init, h, p, chlen) 
  95. {
  96.   //
  97.   // Randwalk uses the transition matrix P, the initial state init,
  98.   // and the number of steps chlen to take a random walk.  The other
  99.   // input variables are required in order to keep the statistics we
  100.   // need to accumulate. 
  101.   //
  102.  
  103.   local (old, new, i, W, total);
  104.  
  105.   old = init;
  106.   W = 1.0;
  107.   total = f[old;1];
  108.   for (i in 1:chlen) 
  109.   {
  110.     new=nextstate(old, P);
  111.     W = W * A[old;new] ./ P[old;new];
  112.     total = total + W * f[new;1];
  113.     old = new;
  114.   }
  115.   return h[init;1] * total ./ p[init;1];
  116. }
  117.  
  118. estimate = function(B, f, h, chlen, niters) 
  119. {
  120.   // Given a matrix equation Bx = f with known coefficient matrix B and
  121.   // known right-hand-side f, the function `estimate' uses a Markov chain
  122.   // technique to approximate the inner product of the known vector h
  123.   // with the unknown vector x (or rather, the inner product of h with the
  124.   // chlen+1st iterate of the fixed point iteration x = Ax+f, where
  125.   // B = I - A).  The input `niters' controls how many random walks are taken.
  126.   // For details of the method, see Chapter 5 of R. Rubinstein's ``Simulation
  127.   // and the Monte Carlo Method.''
  128.  
  129.   local (A, n, i, j, tmp, P, p, total, eta);
  130.   n = B.nr;
  131.   A = eye(n,n) - B;
  132.   
  133.   //
  134.   // Now generate P, the ergodic Markov chain, by scaling the rows of A so
  135.   // that the row-sums are all = 1.  Note: we are making an implicit assumption
  136.   // that the infinity norm of the matrix A is less than one.  (see Rubinstein).
  137.   //
  138.  
  139.   P=zeros(n,n);
  140.   for(i in 1:n) 
  141.   {
  142.     tmp = mysum(abs(A[i;]));
  143.     P[i;] = abs(A[i;])./tmp;
  144.   }
  145.  
  146.   //
  147.   // Generate the initial distribution vector p.
  148.   //
  149.   
  150.   p=abs(h)/mysum(abs(h));
  151.  
  152.   //
  153.   // Now we're ready to take some `walks'.
  154.   //
  155.  
  156.   total = 0;
  157.   for(i in 1:niters) 
  158.   {
  159.     //
  160.     // Generate an initial state for the particle, and then call randwalk.
  161.     // Randwalk will keep track of the statistics we need for each walk.
  162.     //
  163.     j = initstate(p);
  164.     //  printf("Walk no. %d...",i);
  165.     eta = randwalk(A, P, f, j, h, p, chlen);
  166.     //  printf("done\n");
  167.     total = total + eta;
  168.   }
  169.   return total/niters;
  170. }
  171.  
  172. //
  173. // End of file `estimate.r'
  174. //
  175.